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Abstract. Monte Carlo simulation provides a powerful tool for understanding and 
exploring thermodynamic phase equilibria in many-particle interacting systems. Among 
the most physically intuitive simulation methods is Gibbs ensemble Monte Carlo (GEMC), 
which allows direct computation of phase coexistence curves of model fluids by assigning 
each phase to its own simulation cell. When one or both of the phases can be modeled 
virtually via an analytic free energy function [M. Mehta and D. A. Kofke, Mol. Phys. 79, 
39 (1993)], the GEMC method takes on new pedagogical significance as an efficient means 
of analyzing fluctuations and illuminating the statistical foundation of phase behavior in 
finite systems. Here we extend this virtual GEMC method to binary fluid mixtures and 
demonstrate its implementation and instructional value with two applications: (1) a lattice 
model of simple mixtures and polymer blends and (2) a free-volume model of a complex 
mixture of colloids and polymers. We present algorithms for performing Monte Carlo trial 
moves in the virtual Gibbs ensemble, validate the method by computing fluid demixing 
phase diagrams, and analyze the dependence of fluctuations on system size. Our open- 
source simulation programs, coded in the platform-independent Java language, are suitable 
for use in classroom, tutorial, or computational laboratory settings. 
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1. Introduction 

Phase transitions and critical phenomena are topics fundamental to most undergraduate 
and graduate courses in thermodynamics and statistical mechanics [TJ. Aside from their 
intrinsic interest and practical relevance, phase transitions provide a rich conceptual context 
for charting the path from a microscopic Hamiltonian to macroscopic properties of materials 
via partition functions and free energies. Over the past sixty years, molecular simulations 
have provided important insights into the thermodynamic phase behavior of systems ranging 
from simple atomic fluids to complex macromolecular materials. Monte Carlo and molecular 
dynamics methods, in particular, have clarified the subtle interplay between energy and 
entropy in determining stability of competing phases [2J [3]. Monte Carlo simulation has 
been further exploited as an aid in teaching statistical mechanics (U [5j [6]. 

Computer simulations of many-particle systems can readily identify mechanically stable 
(metastable) states, corresponding to local minima in the free energy of the system as 
a function of externally controlled parameters. For fluid systems, external parameters 
may include temperature, pressure, density, and (in the case of mixtures) concentration. 
More difficult is finding the global minimum in the free energy, which is required to 
establish thermodynamic stability. The task is especially challenging near a first-order phase 
transition, where two bulk phases (e.g., vapor and liquid in a one-component system, or 
/L-rich and B-hch phases in a mixture of A and B species) may coexist in equilibrium. 

In simulations of finite model systems, the free energy cost associated with interfaces 
between phases results in hysteresis and tends to inhibit the simultaneous presence of more 
than one phase in a single simulation cell [2]. For this reason, mapping out thermodynamic 
phase diagrams using computer simulation traditionally involves computing the free energy, 
usually via thermodynamic integration of the internal energy, and then performing a 
coexistence analysis via a Maxwell common-tangent construction. An alternative approach, 
which has been applied to both fluid and magnetic systems, is histogram reweighting [6]. 

A more direct and intuitive route to analyzing phase equilibria and calculating densities 
of coexisting phases is the Gibbs ensemble Monte Carlo (GEMC) method. Introduced by 
Panagiotopoulos [TJ El [9j [10], this method models each phase in its own simulation cell. 
By elegantly avoiding both the complication of interfaces and the need for thermodynamic 
integration, the GEMC method provides a computationally efficient and conceptually 
transparent means of computing fluid phase diagrams. The GEMC method has been widely 
applied as a research tool to analyze phase behavior of simple fluids and fluid mixtures [8], 
as well as complex fluids, such as colloid-polymer mixtures [10J [TTJ [12] and charge-stabilized 
colloidal suspensions [13]. To date, however, the pedagogical potential of the GEMC method 
has been relatively less appreciated. 

In this paper, we explore the GEMC method as a tool for introducing students to the 
statistical nature of fluctuations and phase behavior in finite systems. We are inspired by 
a useful variation of the method, proposed by Mehta and Kofke [TJ], which incorporates a 
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thermodynamic model via an analytic free energy or equation of state (EOS). By modeling 
one of the phases explicitly and the other as virtual (via an EOS), the GEMC-EOS method 
gains efficiency over the original GEMC method and is reliable when the EOS of the virtual 
phase is accurately known. 

The GEMC-EOS method was conceived and tested for one-component systems and 
applied to a simple fluid of particles interacting via Lennard- Jones pair potentials. A 
simplified, pedagogically appealing version of the GEMC-EOS method models both phases 
as virtual |14j . The purpose of the present work is to extend this "virtual GEMC method" 
to fluid mixtures and to demonstrate the value of the method in exploring and elucidating 
the statistical nature of demixing. 

The remainder of the paper is organized as follows. In Sec. [21 we briefly review the 
GEMC-EOS method and then describe our extension to fluid mixtures. In Sec. [3j we 
demonstrate the application of the virtual GEMC method to the demixing behavior of two 
model systems: (1) a lattice model of mixtures, and (2) a free- volume model of colloid- 
polymer mixtures. For the latter, we analyze density fluctuations as functions of system 
size and proximity to the critical point. Finally, in Sec. [U we conclude by emphasizing the 
significance of the virtual GEMC method as a classroom and computational laboratory tool. 

2. Methods 

2.1. Gibbs Ensemble Monte Carlo Simulation 

By allowing distinct simulation boxes to exchange particles and volume, but not otherwise 
interact, the Gibbs ensemble Monte Carlo method can efficiently equilibrate coexisting phases 
(see Fig. [p. Each box accommodates one of the phases at equal temperatures, pressures, and 
chemical potentials, with no need for interfaces. Originally introduced by Panagiotopoulos [7] 
to model liquid-vapor coexistence of simple one-component fluids, the GEMC method was 
later extended to mixtures [8]. Since its introduction, the GEMC method has been refined 
and used to map out phase diagrams for a wide variety of systems [9j [TO] . 

Trial moves in Monte Carlo simulations of thermal systems are typically accepted 
with probabilities consistent with the condition of detailed balance [2]. According to this 
condition, the average rate of transition from an old (o) to a new (n) state equals, in 
equilibrium, the average reverse transition rate: 



where P{o) and P(n) are probabilities of finding the system in states o and n and 7r(o — > n) is 
the transition probability between o and n. Assuming that transitions o — > n and n — > o are 
attempted at equal rates, the classic Metropolis algorithm [15] imposes Eq. (TjQ) by accepting 
trial moves with probability [2l [3] 



P(o)ir(o — > n) = P(n)Tc(n — > o) 



(1) 




(2) 
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Figure 1. Schematic illustration of trial moves performed in the Gibbs ensemble Monte 
Carlo method. Top: particle transfers between two boxes in the lattice model of binary 
mixtures (Sec. 13. lj) . Bottom: particle displacements, volume changes, and particle transfers 
in a model of colloid-polymer mixtures (Sec. I3.2[ ). The virtual GEMC method replaces both 
of the boxes by a virtual phase described by a free energy 



In the canonical (constant- NVT) Gibbs ensemble, a microstate of a binary mixture, at 
absolute temperature T, is specified by the positions of all particles (collectively denoted by 
{r}), the volume Vi of each box (i = 1,2), under the constraint of constant total volume 
V = V\ + V2, and the particle numbers Nij of the two species (j = A, B) in each box, with 
constant total numbers N\ = N^i + Nbi and A^ 2 — ^A2 + Nb2- The probability distribution 
for the possible microstates is given by 

PUN*}, {Vi}, T; {r}) cx P id e"^« r » , (3) 

where (3 = l/(ksT), U = U\ + U2 is the internal energy of both boxes, U{ being the internal 
energy of box i, and 



Pid({^},{^})oc 



N Al \N Bl \N A2 \N B2 \ 
is the probability distribution for an ideal mixture. 



(4) 



2.2. Virtual Gibbs Ensemble Monte Carlo: Modeling Phases by a Free Energy 

The GEMC-EOS method replaces one of the simulation boxes by a virtual phase containing 
no explicit particles, but represented instead by a thermodynamic model in the form of a 
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free energy or equation of state. In the virtual GEMC method, both phases are modeled as 
virtual with a prescribed excess free energy F ex , associated with interparticle interactions. 
The probability density for a configuration in which box % has volume Vi and contains Ay 
particles of species j is then given by [cf. Eq. (J3])] 

P oc P id e~ pF ™ . (5) 

According to Eqs. (J2]) and (J^J), the acceptance probability for a trial move from an old 
(o) to a new (n) state is given by the ratio of the corresponding probability densities: 

V(o -)■ n) = min (l, e"^) , (6) 



. Mo, 

where AF ex = F ex (n) — F ex (o) is the associated change in excess free energy. From Eq. (jlj), 
a trial transfer of volume AV from phase 1 to phase 2 (Vi — > V\ — AV, V2 — >■ V2 + AV) is 
accepted with minimum probability 

V m ,= (l-^f] N ' (l + ^f-Y' e-l> AF « . (7) 



In practice, trial moves in ln(Vi/V2) improve efficiency, with acceptance probability [21 E] 

Finally, a trial transfer of a particle, say of species A from phase 1 to 2 (Nai — > Nai — 1, 
A^42 — > ^A2 + 1) is accepted with minimum probability 

trans ^ ^ + j ■ 

In the next section, to illustrate the utility of the virtual GEMC method in exploring 
demixing phenomena, we apply the method to two model mixtures. 



3. Applications of the Virtual GEMC Method to Demixing 

3.1. Lattice Model of Binary Mixtures 

A simple lattice model provides an instructive introduction to phase separation. Consider a 
lattice of N sites, each occupied by a particle of type A or B with volume fractions defined 
as 4>a = Na/N and (pB = Nb/N, the total volume fraction being conserved: 4>a + 4>b = 1 
(see Fig. [1]). The entropy of mixing is given by 

S mix = k B In I N m 1 J - ~h [N A In (f) A + N B In <f> B \ , (10) 

using Stirling's approximation, In A! ~ A In A — A (valid for A ^> 1). Suppose now that 
only nearest-neighbor pairs interact and that correlations between particles can be neglected. 
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Figure 2. Phase diagram for the lattice mixture model. Squares: data from virtual GEMC 
simulations with both phases described by the mixing free energy of Eq. (|13[) . Solid and 
dashed curves: exact binodal and spinodal curves [Eqs. (fT4|) . (fl5]l]. Circle: theoretical 
critical point. Representative tie lines join coexisting phases on the binodal. 



In this mean-field approximation, a particle of type i occupies a given site with probability 
(pi. The internal energy, then independent of configuration, can be expressed as 

U = - (N A <p A e AA + N B <p B e BB ) + N A (j) B e A B , (11) 

where denotes the interaction energy between a pair of particles of species i and j. The 
corresponding internal energy of mixing is 

t/mix = U - - {N A e AA + N B e BB ) = X N A (p B , (12) 

where x = e AB — (^aa + c_bb)/2 is the Flory interaction parameter, which characterizes the 
incompatibility of the two species. Finally, with = <p A , the Helmholtz mixing free energy, 
F mix = [y m i x — TS mix , is approximated (per site) by 

^ = 01n0+(l-0)ln(l-0) + X 0(l-0) . (13) 

The phase behavior of the lattice model predicted by the above mean-field theory is 
easily deduced from the analytic expression for the free energy [Eq. ( TT3] . In equilibrium, two 
phases coexist at equal temperatures, pressures, and chemical potentials, fi A and fi B , of the 
two species. Incompressibility of the lattice and conservation of total volume fraction reduce 
the coexistence criteria to a single condition: \i A = <9-F mix /<90 = 0. The coexistence curve 
(binodal) thus takes the analytic form 

X=-^*(Ar). (14) 
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The inflection points of the free energy define the spinodal curve, 
d 2 F ■ 1 

v 1 mix n 1 /i r\ 

^ = ° X = 20(1^0) ' (15) 

inside of which the mixture is thermodynamically unstable and spontaneously demixes. The 
binodal and spinodal curves terminate at a lower critical point, crit = 0.5, Xcrit = 2, 
below which the mixture is stable. Figure [2] shows the resulting phase diagram, where 
the interaction parameter x plays the role of an inverse temperature. 

In the virtual GEMC method, each phase is governed by the mean-field free energy 
of Eq. f fT3|) . With total volume fraction conserved, the only independent trial moves are 
exchanges of particles between the two boxes. A trial exchange of an A particle in box 1 for 
a B particle in box 2, for example, is accepted with probability 

v u _ „-/3AF mix _ N A1 N B2 -/3AC7 mix nfi x 

(N A2 + l)(N m + l) ' Ubj 

where -F m i x is the total mixing free energy and 

AC/ mix = 2 X (<h - 02 - 1/iVi) (17) 
with (pi = NAi/Ni. For sufficiently large systems, Eq. ( TIBjl can be approximated by 



The Boltzmann factor in Eqs. (fT6|) and (TTST) evidently favors demixing if x > 0, while the 
entropic prefactor always favors mixing. The competition is decided by the magnitude of \- 

As an illustration of the GEMC method with both phases treated virtually, we have 
implemented Eq. ( fTHl) for the lattice model and performed simulations to calculate several 
points on the coexistence curve (binodal). A trial move figuratively flips a coin (i.e., generates 
a random number) to decide, with equal probabilities, whether to exchange an A particle 
for a B particle, or a B for an A, between the two boxes. Choosing N\ = N 2 = 1000, 
and initializing with equal numbers of A and B particles in each box, we equilibrated the 
system for 10 4 MC steps, a step being defined as a trial exchange of every particle. We 
then accumulated data over the next 10 4 steps to calculate average values of 0ai and <pA2 
over a range of x values. Because virtual phases have no explicit particles to displace, these 
simulations are extremely fast, with typical CPU times of a few minutes on a PC, scaling 
linearly with system size. 

Our numerical results, plotted in Fig. |2j agree essentially exactly with the analytic 
expression for the binodal [Eq. (fl4]) ]. thus validating the virtual GEMC method. On 
approaching the critical point, we observed growth of fluctuations and frequent switching of 
phases between the two boxes. Phase switching is easily suppressed, however, by increasing 
the number of particles. In a tutorial setting, students could, for example, probe miscibility 
in different regions of the phase diagram and explore the dependence of fluctuations on 
system size and proximity to the critical point. For this purpose, it is valuable that the 
binodal and spinodal curves are also known analytically. 
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Figure 3. Typical diagnostic data from a virtual GEMC simulation of the lattice mixture 
model [Eq. (|13[) ]: cumulative average volume fraction of A species in each box vs. number 
of Monte Carlo steps during equilibration for interaction parameter x — 2.1. 



To facilitate pedagogical applications on a variety of platforms, we have coded our 
simulations in the Java programming language using the Open Source Physics (OSP) 
library [T8| [T9] . The graphical user interface provided by the AbstractSimulation class 
in the controls package of the OSP library allows the user to easily input parameters and 
start, stop, and step through a simulation. Figure [3] shows typical diagnostic data, collected 
during the equilibration stage, displayed with the OSP frames package. Our code can be 
easily made compatible, however, with any graphics library. 

A closely related model depicts a polymer blend as a mixture of chains of lengths 
(degrees of polymerization) Ma and Mb, whose segments (monomers) fully occupy the sites 
of a lattice. Connecting monomers to form chains reduces their entropic free energy by a 
factor of inverse chain length. Making an analogous mean-field approximation, the Flory- 
Huggins theory for this model yields p2] 

^ = j£w + l^Mi-*) + x#(i-*), < 19 > 

where iV^ now represents the number of A chains (rather than A segments) and = 
Na_Ma/N. In the symmetric case (Ma = Mb = M), the phase diagram is identical to 
Fig. [21 but with \ replaced by M\ on the vertical axis. Asymmetric blends (Ma ^ Mb) 
display much richer phase behavior. The Flory-Huggins mixing free energy [Eq. (I19p ] could 
be used to explore, via the virtual GEMC method, demixing of polymer blends. 
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3.2. Colloid- Polymer Mixtures 

As a second application of the virtual GEMC method, we consider a widely studied system 
from soft matter physics - a suspension of colloidal particles mixed with free (nonadsorbing) 
polymers [201 122]. The classic Asakura-Oosawa-Vrij (AOV) model [231 121] idealizes the 
colloidal particles as hard spheres, monodisperse in radius R c , and the polymer coils as 
effective spheres with a radius R p equal to the average radius of gyration. The polymers are 
modeled as mutually ideal (noninteracting) , but impenetrable to the colloids, with which 
they have hard-sphere interactions. Although real polymer coils fluctuate in size |25j, the 
AOV model portrays the effective polymer spheres as monodisperse in size (see Fig. [T)). The 
size ratio q = R p j R c is then the one model parameter that distinguishes different mixtures. 

The thermodynamic state of the system is specified by the total volume V and numbers 
of colloids and polymers, N c and N p , with respective number densities p c = N c /V and 
P P = Np/V and volume fractions C = (47r/3) p c R^. and <p p = (4ii/3)p p R p . In the Gibbs 
ensemble, the particle numbers in the two boxes are denoted N c i, N c2 , N p \, N p2) and 
constrained by N C \+N C 2 = N c and N P \+N P 2 = N p . With only hard interparticle interactions, 
the thermodynamic state is independent of temperature, there being no energy scale. For 
contact with experiments, it is helpful to imagine the system in osmotic equilibrium with 
a reservoir of pure polymer, which exchanges polymers with the system to maintain fixed 
polymer chemical potential. The reservoir density p p plays the role of an inverse temperature. 

To describe the phase behavior of the AOV model of colloid-polymer mixtures, 
Lekkerkerker et al [26] have developed a free-volume theory by expressing the Helmholtz 
free energy density (to within a constant) as 

f{<t>c <t> P ) = k B Tp c (In (j) c - 1) + f hs {(j>c) + fp(<frc <Pp) ■ (20) 

The first two terms on the right side are the colloid ideal-gas and hard-sphere excess 
free energy densities, the latter being accurately approximated by the Carnahan-Starling 
relation [28J: 

C (4 - 

If colloid-polymer correlations are neglected (mean-field approximation), the polymer free 
energy density can be approximated by that of an ideal gas of polymers confined to the free 
volume, i.e., the volume not excluded by the hard-sphere colloids: 



pfhsWc) = Pc— T\2~ ■ ( 21 ) 



fp{<Pc,<Pp) = k B Tp p 



In 



(22) 



Here a(<f) c ) is the free-volume fraction of polymers amidst colloids, which is reasonably 
approximated by scaled-particle theory: 

3 



a{<p c ) = — exp 

1 — Or 



J2 C ^ m ) > ( 23 ) 



m=l 
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Figure 4. Phase diagram for the AOV model of colloid-polymer mixtures with size ratio 
(polymer radius of gyration over colloid radius) 3=1. Squares: data from virtual GEMC 
simulations with both phases described by free energy from free-volume theory. Solid and 
dashed curves: binodal and spinodal curves from Maxwell common-tangent construction. 
Circle: theoretical critical point. Tie lines join coexisting phases on the binodal. 

where 7 = C /(1 - (j) c ),C x = 3q + 3q 2 + q 3 ,C 2 = (9q 2 /2) + 3q 3 , and C 3 = 3q 3 . For 
ideal polymers, equality of polymer chemical potentials in the system and reservoir implies 
4> p = 0ptt(0 c ), where (p r p = (4n/3)p 7 p R 3 is the volume fraction of the polymer reservoir. 

Given the analytic expression for the free energy [Eqs. ( !20|) -(l23l)]. it is straightforward 
to explore phase stability of the AOV model predicted by the mean-field free- volume theory. 
The demixing binodal can be calculated from a Maxwell construction that equates pressures 
and chemical potentials of coexisting phases, while the spinodal is defined by the inflection 
points of the free energy. Figure H] shows the resulting fluid demixing phase diagram for 
mixtures with a size ratio q = 1, for which the polymer coils are impenetrable to the colloids. 
The colloid-rich and colloid-poor branches of the binodal correspond, respectively, to colloidal 
"liquid" and "vapor" phases. 

For simulations, the free energy expressed by Eqs. ( I20l) -()23l) can be recast in the form 

/(0c, <t>p) = /id(0c, 0p) + /cxOc <t>p) , (24) 

comprising an ideal-gas free energy density, 

/3f id ((f) c ,(f) p ) = p c (\n(p c - I) + p p (In (p p -l) , (25) 
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and an excess free energy density 

PfMc P ) = Pfhs(<Pc) ~ P P lna(0 c ) . (26) 

The acceptance probabilities for the trial moves in the virtual Gibbs ensemble are given by 
Eqs. (jHJ) and (J9]) combined with Eq. ( 126]) . A transfer of volume AV from phase 1 to 2 is 
accepted with probability 

AV\ Nl+1 { AV^ N2+1 



[^^) e"^-, (27) 

while transferring a colloid from phase 1 to 2 is accepted with probability, 

Vtrmis = Yl J^£L_ . (28) 

Vi N c2 + 1 y 1 

The last factor on the right side of Eqs. ( |27|) and ( 128|) may be expressed as 
-/?AF ox _ f ai(n)\ Npl ( a2{n)\ Np2 « AFh 



To illustrate our implementation of the virtual GEMC method for the AOV model, 
we have performed simulations of mixtures with a size ratio of q — 1. For efficiency, we 
pre-computed fh s and a vs. C and stored these data in a lookup table. Initializing the 
system with equal numbers of particles and equal volumes in the two boxes, we fixed the 
total number of colloids at N c = 2000, chose the volume to give an average colloid volume 
fraction of <j) c = 0.1, and adjusted the polymer volume fraction by varying the total number 
of polymers over a range 2500 < N p < 10 4 . After equilibrating for 10 4 steps, we accumulated 
statistics over the next 10 4 steps to calculate average volume fractions of each species in each 
box. A step is here defined as one trial volume exchange and a trial transfer of every particle. 
In attempting a particle transfer, we first randomly chose a box (1 or 2), then randomly chose 
a species of particle to be transferred to the other box. 

The results of our simulations are plotted on the demixing phase diagram in Fig. HJ 
Each run generated a pair of points on the binodal - one on the liquid branch and one on 
the vapor branch. Theory and simulation again agree very closely, further validating the 
method. As with the lattice model, students could, in a computational laboratory setting, 
explore phase stability in different parts of the phase diagram and explore the variation of 
fluctuations with system size and proximity to the critical point. Although the free- volume 
theory of the AOV model does not yield the binodal and spinodal curves in analytic form, the 
calculation of these curves by a Maxwell construction on the free energy is straightforward. 

Figure [5] shows a typical diagnostic trace near the beginning of a run in the unstable 
region, revealing significant fluctuations in the densities of the coexisting phases. To analyze 
these variations, it is helpful to define the root-mean-square (rms) fluctuations 
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Figure 5. Typical diagnostic data from virtual GEMC simulation of AOV model of colloid- 
polymer mixtures: cumulative average volume fractions of colloids in each box vs. Monte 
Carlo steps during equilibration for size ratio q = 1 and total volume fractions <j) c = 0.1, 
and 4> p = 0.3. 

where cj>i is the mean volume fraction of species i (i = c, p) and cj^p is the volume fraction of 
species i sampled at Monte Carlo step j. Figure [6] shows the relative rms fluctuations, cr c /(fi c 
and a p /(f)p, for two sets of systems - one relatively small, with colloid number N c = 100 and 
total polymer numbers in the range N p = 6 x 10 2 -10 3 , and the other 10 times larger, with 
N c = 1000 and N p = 6 x 10 3 -10 4 - computed from runs of 2 x 10 5 MC steps. To minimize 
correlations between successive samples of the volume fraction, we spaced the samples by 
intervals of 1000 steps. 

As seen in panel (a) of Fig. El fluctuations grow upon approaching the critical point 
and as the system size is decreased. Panel (b) replots the data from panel (a), but with 
fluctuations in the larger systems now scaled by a factor of a/10 (square-root of system 
size ratio). The collapse of the two data sets upon scaling demonstrates the well-known 
result that relative fluctuations scale with the inverse-square-root of the particle number: 
o~i/(j>i ~ 1/yN. As an exercise, students could test this scaling property by performing 
simulations over a range of system size. 
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Figure 6. Root-mean-square (rms) fluctuations in densities of coexisting phases, 
corresponding to the phase diagram in Fig. |4] (a) Relative fluctuations a c /4> c and a p j 1 <p p 
for systems with polymer numbers in the range N p = 6 x 10 2 -10 3 (squares) and for systems 



with 10 times as many particles, N p — 6 x 10 3 -10 4 (circles), 
larger-system fluctuations now scaled by a factor of y/lO. 



(b) Same data, but with 
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In summary, we have extended to fluid mixtures a version of the Gibbs ensemble Monte Carlo 
simulation method that models both phases virtually via a thermodynamic equation of state 
or free energy. As illustrations, we have applied the method to demixing in two models for 
which, in a mean-field approximation, analytic free energy expressions are known: a lattice 
model of simple mixtures and a free-volume model of colloid-polymer mixtures. For both 
models, we have validated the method by computing fluid demixing phase diagrams that 
closely agree with those calculated from Maxwell common-tangent constructions. 

As a computational method, virtual GEMC has two main virtues. First, for finite 
systems whose equation of state may be known approximately from experiment, but 
whose interparticle interactions remain unknown, the method provides an alternative to a 
Maxwell construction that incorporates the impact of fluctuations on demixing. Second, 
and perhaps more importantly, the virtual GEMC method has pedagogical value as a 
tool for demonstrating the statistical nature of phase behavior and for allowing rapid 
exploration and analysis of fluctuations in finite systems. Our Java simulation programs 
can be readily adapted for use as classroom demonstrations or as numerical "experiments" 
in a computational laboratory. Finally, it should be straightforward to generalize the virtual 
GEMC method to other models and to mult i- component mixtures. Future applications, for 
example, may explore demixing phase diagrams of ternary oil-water-surfactant mixtures. 
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